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Abstract. We present numerical models of hydromagnetic instabilities under the conditions prevailing in a stably stratified, 
non-convective stellar interior, and compare them with previous results of analytic work on instabilities in purely toroidal fields. 
We confirm that an m = 1 mode (‘kink’) is the dominant instability in a toroidal field in which the field strength is proportional 
to distance from the axis, such as the field formed by the winding up of a weak field by differential rotation. We measure 
the growth rate of the instability as a function of field strength and rotation rate ff, and investigate the effects of a stabilising 
thermal stratification as well as magnetic and thermal diffusion on the stability. Where comparison is computationally feasible, 
the results agree with analytic predictions. 
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1. Introduction 


Magnetic fields probably play a significant role in the internal 
rotation of stars. Even a relatively weak magnetic field is suf¬ 
ficient to couple different parts of the star and maintain a state 
of nearly uniform rotation. For the interior of the present Sun, 
for example, a field of less than 1 gauss would be able to trans¬ 
mit the torque exerted by the solar wind through the interior 
dMestel 1953K The observed rotation in the core of the Sun 
( jChaplin et al. 200T| is quite uniform, suggesting that a mag¬ 
netic field of this order or larger may actually be present. The 
progenitors of white dwarfs and supernovae go through giant 
stages in which the envelope rotates very slowly. The degree 
of coupling between core and envelope by a magnetic field in 
this stage will determine whether the rotation rates of pulsars 
and white dwarfs are just a remnant of the initial rotation of 
their progenitors, or if a secondary process must be responsible 
( jSpruit & Phinney 1998l|Spruit 1998} . 


Models of gamma-ray bursts in which the central en¬ 
gine derives from the rapidly spinning core of a massive star 
( |Woosley 1993||Paczynski 19 981 also depend on the ability of 
the core to keep its high angular momentum for a sufficient pe¬ 
riod of time, in the face of magnetic spindown torques exerted 
by the slowly rotating envelope (|Heger et al. 2000||. 


The uniform rotation of the solar core may be due to a mag¬ 
netic field, but this field’s origin, configuration and strength are 
not known. By analogy with the magnetic A-stars, one might 
speculate that a ‘fossil’ magnetic field could exist in the core 
of the Sun i jCowling 1945) [Braithwaite & Spruit 2004 1 . Since 
no net field is seen at the surface (averaged over the solar cy¬ 
cle), the radial component of such a fossil would however have 
to be weak - of the order of a gauss or less. A field weaker 


than about this 1 G will quickly wrap up into a predominantly 
toroidal field, under the action of the remaining (weak) differ¬ 
ential rotation in the core. 

The predominantly toroidal magnetic field resulting from 
this process will not increase in strength arbitrarily. Eventually, 
the energy density in the field will become large enough that a 
magnetic instability will set in. 

Analytic work, (e.g. [Tayler 1973| , shows that any purely 
toroidal held should be unstable to instabilities on the mag¬ 
netic axis of the star (pinch-type instabilities, under the in¬ 
fluence of the strongly stabilising stratification in a radiative 
stellar interior, or ‘Tayler instabilities’ hereafter). The growth 
rates of these instabilities are expected to be of the order of 
the time taken for an Alfven wave to travel around the star on a 
toroidal held line. This is very short compared to the evolution¬ 
ary timescale of the star. In a star like the Sun, for example, with 
a held of 1000 gauss, the growth timescale ry'ATrp/B would be 
of the order of years, if r = Rq/2 is taken, and p = 1.3g/cm 3 . 

A magnetic held of this type can also be subject to other 
instabilities, such as the magnetic buoyancy ( [Parker 195~5l 
IGilman 19701 lAcheson 19781 and magnetic shear instabilities 
d Velikhov 19591 lAcheson 1971 |Balbus & Hawley 1992) . As 
was shown by Spruit (1999), the Tayler instability will be the 
first to appear as the strength of the toroidal field is increased. 
This is because with the magnetic buoyancy instability, as with 
all instabilities where displacements in the vertical are neces¬ 
sary, the stratification provides a strong stabilising force. The 
same is the case with the magnetic shear instability, whose ef¬ 
fect in a stellar interior is very limited in comparison to its ef¬ 
fect in accretion discs. In contrast, the Tayler instability occurs 
on the magnetic axis, where the magnetic field is perpendicular 
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to gravity and the displacements caused by the instability are 
also perpendicular to gravity. 

In this paper we aim to test numerically the instability 
mechanism, and to verify that the predictions of the analytic 
work are relevant: that they cover all instabilities actually 
present in a system consisting of a predominantly toroidal field 
in a stable stratification. Much of the analytical stability anal¬ 
yses have been done under a local approximation. This can be 
shown to be exact for the case of adiabatic instabilities in a non¬ 
rotating star, but not for the more interesting cases in which 
rotation and the effects of magnetic and thermal diffusion are 
taken into account. Though it is not expected that major insta¬ 
bilities have been missed, numerical simulations can provide 
an important check. 

It is much less certain how the magnetic field evolves once 
instability has set in. In a scenario developed by Spruit (2002), 
it is argued that the instability will lead to self-sustained dy¬ 
namo action. The field remains predominantly toroidal, subject 
to decay by Tayler instability, but is continuously regenerated 
by the winding-up of irregularities produced by the instability. 
This scenario has been applied in stellar evolution calculations 
of the internal rotation of massive stars by Heger et al. (2003) 
and Maeder and Meynet (2003). 

The balance between wrapping-up by differential rotation 
on the one hand and the destruction of the toroidal field by 
Tayler instabilities on the other determines the strength and 
configuration of the field, and the the rate at which it transports 
angular momentum through the star. 

The long-term goal is therefore to investigate the non-linear 
development of the instability. With 3-D numerical simulations 
we can determine how quickly an initial toroidal field decays 
(by reconnection across the magnetic axis), and to determine 
the type of magnetic field that is maintained by differential ro¬ 
tation in a stably stratified star, under the action of magnetic in¬ 
stabilities, and to develop from this a quantitative theory for the 
angular momentum transport by magnetic fields in stars. This is 
beyond the scope of this paper, but is looked at by Braithwaite 
(2005), where the operation of this differential-rotation driven 
magnetic dynamo is demonstrated. 

2. The nature of the instabilities 

We first consider adiabatic instabilities, that is, ignoring the ef¬ 
fects of viscosity and of thermal and magnetic diffusion. The 
instabilities of a magnetic field in a stable stratification then 
depend on three parameters: the field strength, some measure 
of the stability of the stratification, and the rotation rate of the 
star. For stars rotating well below their critical (maximal) rate, 
and for the expected relatively low field strengths, the relative 
strength of the parameters is expressed by the ordering 

TV 3> ft and N ft a, (1) 

where N is the buoyancy frequency, ft is the rotational ve¬ 
locity, and ft a is the Alfven frequency given by va/R = 
B/(Ry/4np), R is the radius of the star and va is the Alfven 
speed. This is the same, in effect, as saying that the thermal 
energy density is much greater than both the rotational and the 
magnetic energy densities. 




equilibrium m = 0 


m = 1 


m = 
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Fig. 1. The physical form of the instability, for modes m = 0, 
m = 1 and m = 2. Above, co-moving surfaces are drawn, with 
the magnetic held represented by the lines with arrows. The z 
and vj axes are marked on the diagram on the left; gravity (if 
present) acts along the z axis. Below, cross-sections of these 
co-moving surfaces are plotted, with the dotted circles repre¬ 
senting equilibrium. 
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In cylindrical coordinates (ru, <l>. z ) where the 2 axis is par¬ 
allel to the magnetic axis, the displacements from equilibrium 
are of the form 


S 


(3) 


The shape of the unstable displacements is shown in Fig. 
13 for modes with m = 0,1, 2. The to = 1 mode is known 
as the ‘kink’ instability, and is seen in twisted solar coronal 
loops and twisted rubber bands. This instability sets in once 
the magnetic held component along the flux tube (vertical in 
the geometry used here) falls below some critical value in rela¬ 
tion to the toroidal component. In our case of a predominantly 
toroidal held, resulting from winding up by differential rota¬ 
tion, the vertical component is neglected: it is a limiting case. 

It is predicted that for an instability to occur, we must have 
( |Tayler 1957) 

m 2 

p > —-1 (m/0) and p > 1 (ra = 0). (4) 
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It is now useful to consider the likely value of p in a 
differentially-rotating star. The toroidal field wound up by dif¬ 
ferential rotation will be: 


dB ^ 
dt 


p> init ^ , pinit 

= + "to’ 


(5) 


where By 1 " and B'"" are the initial field components. We ex¬ 
pect both spatial derivatives of f l to tend to some finite value 
on the z-axis, whether the differential rotation is being driven 
by braking on the surface or by evolutionary changes in the 
moment of inertia, etc. There is no obvious reason why £?‘ nit 
should go to zero on the axis, so the first term on the right-hand 
side of Eq. □will be proportional to vj. The behaviour of B'Jl" 
near the z-axis is less certain, and it may have some dependence 
on vo so that the second term goes as vo 2 or higher, but overall 
the toroidal field component B & will be proportional to vo, i.e. 
P= 1- 

So if we use p = 1, the only unstable mode is m = 1; 
to = 0 and to = 2 are marginally stable. If p = 2, then we 
should expect the to = 0, to = 1 and m = 2 modes all to be 
unstable. 

Only high values of the vertical wavenumber n are unsta¬ 
ble. This is due to the dominant effect of the stable stratifi¬ 
cation. In order to minimise the energy lost against the stable 
buoyancy force, the vertical displacements have to be small 
compared to the horizontal displacement. Since the displace¬ 
ments also have to be nearly incompressive (otherwise energy 
is lost in doing work against the pressure force) this implies 
that the vertical length scale is small. 

If we define a local Alfven frequency 

oja = v A /vo = B j(w\J 4:irp), ( 6 ) 


then we expect the growth rate for a non-rotating star to be 
( [Tayler 1973| and fGoossens et al. 19811 

er ~ wa (fi < wa). (7) 


expect stability at all values of m. [We note that this is a pecu¬ 
liarity of the adiabatic case: when magnetic and thermal diffu¬ 
sion are included, the conditions on p for instability reverts to 
0 (cf. Appendix in Spruit 1999 and below).] 


2.1. The effect of magnetic diffusion 

As mentioned, there is a lower limit on the vertical wavenum¬ 
ber of an unstable mode, caused by the work which must be 
done against gravity to move gas in the vertical direction. There 
is also an upper limit to n, caused by diffusion: the unstable per¬ 
turbations diffuse away at a rate of order r/n 2 , so an instability 
whose intrinsic (adiabatic) growth rate is less than its decay rate 
by diffusion will be smothered. The range of unstable radial 
wavenumbers for the m = 1 mode is therefore, approximately, 
JAcheson 1978l|Spruit 1999| >: 


<7 2 

- > n z > 

1 


N 2 


, 2 r 2 ‘ 


(9) 


where r is some measure of the size of the field in the hori¬ 
zontal direction. This is the condition in the absence of thermal 
diffusion (re = 0). A complication arises here because it turns 
out that in the presence of both thermal and magnetic diffusion, 
the limit k/t] —> 0 is a singular one with respect to the insta¬ 
bility condition: there exist double diffusive instabilities which 
are absent in the case re = 0. Since there is always some ther¬ 
mal diffusion due to numerical effects, we should not expect 
the simulations to reproduce 0 accurately. This is discussed 
further below. 

For a p = 1 field with B = (0, Bqvo/voo , 0) in the slowly- 
rotating case Q <C u>a we have unstable wavenumbers n, and 
hence instability, if (using Eq. 0 ) 


wa > 



( 10 ) 


The instability condition is local (in a meridional plane): if 
it is satisfied at some point (vo, z), an unstable eigenfunction 
can be fit into a small (for adiabatic instability: infinitesimal) 
region around this point (Tayler 1973). The instability is thus 
of an ‘interchange’ type. This property holds by virtue of the 
assumption of a purely toroidal field. It greatly simplifies the 
analytical stability analysis, and allows a detailed treatment of 
the effects of diffusion and viscosity dAcheson 19781) . In the 
azimuthal direction, the instability is global, since only low az¬ 
imuthal orders to are unstable. Connected with this is the fact 
that the typical instability time scale (if conditions for instabil¬ 
ity are satisfied) are always of the order of the time it takes an 
Alfven wave to travel around the star in the azimuthal direction. 

If the star is rotating with the rotation axis parallel to the 
magnetic axis, a stabilising effect is produced so that instead of 
0 we have ( Pitts & Tayler 1986$ 

ui 2 

p > + 1 (to ^ 0) and p > 1 (to = 0). (8) 

in the limiting case where the rotational velocity is much 
greater than the magnetic. In the particular case of p = 1, we 


Instability is thus suppressed at low field strengths. The in¬ 
stability condition is a function of the meridional coordinates 

(w,z). 

2.2. The effect of thermal diffusion 

Thermal diffusion has the effect of reducing the stabilising ef¬ 
fect of the density stratification, allowing instability for a larger 
range of vertical wavenumbers. This is a general effect, not 
only for the magnetic instabilities discussed here. It was first 
noted in the astrophysical context by Zahn (1974). It is im¬ 
portant in stars since the thermal diffusion (measured by the 
diffusivity, with units area per unit time) is much faster than 
other damping effects (like magnetic and viscous diffusion). In 
this case there exist length scales large enough that these other 
damping effects are negligible on the time scale on which the 
instability operates, but at the same time small enough that ther¬ 
mal diffusion can wipe out the temperature fluctuations due to 
vertical displacement against the stable thermal stratification. 
This led to Zahn’s (1974,1983) formulation of shear instability 
in stellar interiors, widely used in stellar evolution calculations. 
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The effect can be incorporated by replacing the buoyancy 
frequency frequency TV by a value N given by 


N 2 


N 2 

1 + tj/tt ’ 


( 11 ) 


where tj is the timescale of the instability, i.e. 1/a, and tt is 
the thermal timescale, equal to 1 /m 2 . [This is correct by order 
of magnitude for all instabilities, and can be made more exact 
for any given one]. This should manifest itself in a decrease 
of the minimum unstable wavenumber given in (|9j above; it 
has exactly the same effect as a reduction in the acceleration of 
gravity g. 


In the calculations, we have left out the ohmic heating term 
77 J 2 /p from the energy equation d 1 5b . In a real star, there are 
other sources and sinks of energy (radiation) which together 
determine its thermal equilibrium structure. Since the magnetic 
energy B 2 is small, thermal changes due to Ohmic dissipation 
are slow compared with the (Alfvenic) time scale of interest 
here, it is more consistent to leave out the gradual heating by 
Ohmic diffusion together with these other sources of slow ther¬ 
mal change. 


3. 3D MHD simulations 


2.3. The model 

To simulate conditions appropriate to a stratified stellar inte¬ 
rior, without simulating the star as a whole, we make use of the 
fact that the Tayler instability is always present in at least some 
region near the magnetic axis. It has a low threshold, being sup¬ 
pressed only at very low field strengths where the combined ef¬ 
fects of vertical stratification and magnetic diffusion suppress 
the instability (cf. Eq. illOH . We use a local plane-parallel ap¬ 
proximation of the star, with the direction of gravity parallel to 
the magnetic axis. This is a good approximation for a small re¬ 
gion near the axis at the low field strengths of primary interest, 
since the (meridional components of the) unstable wavenum¬ 
bers are then high (cf. Eqs. and 03 )- 

The ideal gas equation of state, measuring temperature in 
units such that the molar gas constant divided by the molecular 
mass is unity, is 

P = pT and e =-^—. (12) 

7-1 

The momentum equation: 

—- = —-VP + -J xB + g + 2uxfi. (13) 

D tp p 

Conservation of mass: 

§? = -' Vu - (14) 

Specific internal energy: (thermal conduction and Ohmic heat¬ 
ing) 

De 1 , 1 9 

— = — TV.u + r yK-V.(pVe) + p—J 2 . (15) 

ut p p 

Electric current: 

J = V x B, (16) 

and the induction equation: 

— = -Vx(t]J-uxB), (17) 

at 

where in the above equations, the velocity field is denoted 
by u, the differential operator D/Dt = d/dt + u.V, 7 and 
k are the magnetic and thermal diffusivities respectively, e is 
specific internal energy and g is gravitational acceleration. We 
ignore fluid viscosity. 


3.1. The numerical code 

We use a three-dimensional MHD code developed by Nordlund 
& Galsgaard (1995). The code uses a Cartesian coordinate sys¬ 
tem. This has the advantage, over a code using cylindrical co¬ 
ordinates, not only that the code is significantly simpler, but 
that it avoids the coordinate singularity on the axis, known to 
cause serious problems in essentially all grid-based methods in 
cylindrical or spherical coordinates. The disadvantage is that 
an intrinsically round peg (the toroidal field) has to be fitted 
into a square hole. The penalties are some waste of computing 
time (the grid points in the corners being unused), and a small 
amount of startup noise because a good initial equilibrium state 
is harder to construct. 

The code uses a staggered mesh, so that variables are de¬ 
fined at different points in the gridbox. For example, p is de¬ 
fined in the centre of each box, but u x in the centre of the 
face perpendicular to the x-axis, so that the value of x is lower 
by ^dx. Interpolations and spatial derivatives are calculated 
to fifth and sixth order respectively. The third order predictor- 
corrector time-stepping procedure of Hyman (1979) is used. 

For numerical stability, the code contains high-order damp¬ 
ing terms. These were however switched off in this study - we 
were interested in the early (and therefore linear) stage of in¬ 
stability growth, before displacements and velocities become 
large enough to necessitate any stabilisation. 

3.2. Definition of the stratification 

We want to investigate, amongst other things, which vertical 
wavenumbers are unstable. This is achieved most easily if the 
growth rates of the instabilities, and indeed all timescales, are 
independent of z. To this end T has to be made independent of 
z, so that the sound crossing time over the width of the box is 
independent of height. The stratification is thus approximated 
by that of an isothermal atmosphere. This is sufficient for the 
present purpose, since it allows the inclusion of the essence of 
the stabilising effect of buoyancy. 

The equilibrium magnetic field B is then made dependent 
on z in such a way that the ratio of thermal and magnetic energy 
densities, 3 = ep8n/B 2 , and therefore the Alfven speed, is 
independent of z. 
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3.3. Initial magnetic configuration and equilibrium 

Let (zu, (j>) denote polar coordinates on a horizontal plane. We 
consider a purely azimuthal field of the form 


B = Bn I 


(18) 


where wq is a length constant, set to a quarter of the horizon¬ 
tal size of the computational box, and B 0 is a constant which 
determines the strength of the magnetic field. The exponen¬ 
tial containing z keeps the value of 3 independent of z (see 
Sect. 13.4.21 . The instability is in the centre, where w/wq <C 1. 
The exponential containing vj is present to make the axisym- 
metric field fit inside a square computational box. As the field 
evolves during the instability, it will spread outwards so that 
some margin has to be left empty around the region contain¬ 
ing the field. The results show that this has been successfully 
avoided, at least during the linear phase of growth. 

To ensure that the initial state is close enough to equilib¬ 
rium to convincingly measure the growth of the instability, the 
pressure, density and temperature fields have to be adjusted to 
the magnetic force brought about the field given in (II 81 , 

Let 

p(vj, z) = poW(vj)e~ 9Z ^ T ° and T(w) = T 0 D(w). (19) 

The equilibrium conditions (Eqs. (II 21 - o with dt = v = 
rj = 0 ) are then satisfied for the axisymmetric field given by 
d provided that 


W{w) = 1 + 
D{vj) = 1 — 


B 0 2 ct 2 /ZUn 

4po?o 

2tJ7 2 /tx7g 

1 + Wo 

B o 


and 


for the case p = 1 , or, if p = 2 , 


W(m) = 1 + 


D(vj) = 1 — 


Bq 

SpoTo 


1 + 2- 77 


Aw a /rojj 


- 2 vj 2 /-c 


1 + 2 1 


8poTp 2zcr 2 /-CC7? 
B‘ 2 e 


( 20 ) 

( 21 ) 

and ( 22 ) 

(23) 


Magnetic diffusion will destroy the equilibrium even in the 
absence of dynamic instability. This should not be problematic 
as long as the timescale over which the instability becomes vis¬ 
ible is shorter than the timescale over which the equilibrium is 
lost in this way. This condition is comfortably fulfilled in the 
simulations presented here, because the instability operates on 
a small length scale. Magnetic diffusivity of a value which is 
relevant for the instability has little effect on the much larger 
length scale of this initial field configuration. 


3.4. Boundary conditions 

Under the conditions prevailing in stellar interiors, the expected 
length scale of the instabilities is very small (at least in the ra¬ 
dial direction). As in the case of hydrodynamic instabilities, 


therefore, different parts of the star are in effect disconnected 
from each other as far as the instability is concerned. The com¬ 
putational box can then be taken to cover only a small part of 
the star, as is also done in studies of magnetic instability in thin 
accretion discs (e.g. Hawley et al. 1995). This raises a prob¬ 
lem with the boundary conditions, since the box boundaries 
are special locations in the computational volume that have 
no counterpart in the real star. With rigid boundaries, spurious 
phenomena such as boundary layers would affect the results. 
Such effects are minimised by using periodic boundary condi¬ 
tions: there are then no special locations in the computational 
box. 


3.4.1. Horizontal 

Periodic conditions in the x and y directions are implemented 
by copying unknown values of P, T, p and v outside the box 
boundaries from their values at x ± L, y ± L where L is the 
width of the box. For the magnetic field, the sign of B is re¬ 
versed in the copying process, so that there are no current sheets 
at the boundaries. In any case, with the field used, the values at 
the boundaries are very small compared to the values at the 
centre of the box. 


3.4.2. Vertical 

In the vertical direction periodic conditions are also possible 
by making use of an invariance of the equations for the case of 
an isothermal stratification. By this symmetry, a shift in height 
is equivalent to multiplication of the variables by constant fac¬ 
tors. In this way the variables can be scaled from the bottom to 
the top of the computational box. Thus, when interpolating or 
differentiating across the top and bottom of the computational 
box, we multiply or divide the density by a factor exp(L z /H) 
(where L z is the height of the box and H is the scale height). 
These ‘scaled periodic’ boundary conditions are the appropri¬ 
ate equivalent of periodic conditions for a stratified medium. 

Unfortunately, it turns out that this procedure does not work 
in its simplest form. Small perturbations launched by the ini¬ 
tial disequilibrium propagate vertically, as a mixture of sound 
waves and internal gravity waves. The upward propagating 
components of a sound wave will grow (in the limit where the 
wavelength is much smaller than the scale height) at a rate of 
the order of the buoyancy frequency (differing from it by a fac¬ 
tor 1.02 if 7 = 5/3). If the wave leaving the top of the box 
is fed back in at the bottom, this amplification goes on indefi¬ 
nitely. 

Though the initial perturbations can be made small, the ex¬ 
ponential growth of the upward-travelling waves causes a prob¬ 
lem if this growth is faster than that of the magnetic instabil¬ 
ity under study. Since we want to operate in the regime where 
N wa. this will indeed be the case. Engineering fixes like 
artificial damping terms designed to affect sound waves would 
not greatly help, since internal gravity waves will also grow in 
this way at a comparable rate. 

Many methods were tried to extinguish these waves. The 
most satisfactory solution found was to divide the volume ver- 
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tically into two halves, with gravity pointing downwards and 
upwards in the top and bottom halves of the box respectively. 
A sound wave is still amplified within one half, but will not 
see any net increase in amplitude over the course of a journey 
through the whole box. The physics of the magnetic instability 
is replicated in the two halves, so that the penalty is a factor of 
two in computational cost. 


4. Results 

To test the analytic predictions and the behaviour of the code a 
series of test cases was studied, starting with a simple unstrat¬ 
ified medium with adiabatic conditions. Physical ingredients 
are then added concluding with the case of most interest for 
a stellar interior, in which stratification, magnetic and thermal 
diffusion are included. 


3.5. Initial perturbations 

All of the runs were started with a small perturbation to the ve¬ 
locity field. This could either be given at all values of m and n, 
or just to one particular mode. There was some uncertainty as 
to the precise form of the perturbation required; how it should 
vary as a function of w, whether it should include a vertical 
component and in what sense the phases of the two horizontal 
components should vary with z were all unanswered questions. 
Several different formulations were tried, producing, reassur¬ 
ingly, virtually identical results. For the runs appearing here, 
the initial perturbation had its maximum at vj = 0 and died 
away at increasing w so that at vj = wq it was almost zero, the 
vertical component of the velocity perturbation was zero and 
for the to A 0 modes, the phase in the <b direction went from 
0 to 27r over the course of one vertical wavelength so that the 
cross-section remained the same at all z but was rotated about 
the axis once per wavelength. The initial perturbation is there¬ 
fore of the form: 


v 


E 

m=0,n 


V? 


VJ 

— exp 
VjQ 



cos (nz)e ZJ 


+ 


E V ™ eX P 

m^0,n 




[cos(nz — m,</)e CT 


+ sin(nz — m(j>)e^]. (24) 


The amplitude of the initial velocity perturbation V/J was 
given a value 3 x 10" 5 , which corresponds to a fraction 2 x 
10 ~ 5 of the sound speed or 2 x 10 -4 of the Alfven speed. This 
value is large enough that the numerical perturbation is small 
in comparison, while being small enough to follow the linear 
phase of instability growth for several growth time-scales. 


3.6. Accessible parameter values 

The separation of time scales wa -C N and (> <7 A' appro¬ 
priate for stars is also convenient for the analytical stability 
analysis. Such widely differing time scales are problematic for 
numerical simulations, however, since the time step will be set 
by the shortest time scale, the sound crossing time (which is of 
the order of 1/N), while the phenomena of interest happen on 
the slower magnetic time scale. In the first set of calculations 
reported below we have set Q = 0, and the ratio oj a /N is of 
order 0.01 — 0.1. For such values, the predicted growth rates 
are already close to their asymptotic (ojx •C N) value, so not 
much would be gained from more expensive calculations with 
lower field strengths. The effect of rotation is investigated in 
Sect. IQ 


4.1. Dependence of stability on p 

Two initial field configurations, p = 1 and p = 2, were tried. 
A resolution of 48 x 48 x 48 was used for these runs, and the 
computational box is a cube of side 27 t. A zero value was used 
for g , and Bq = 1.28. As in all runs 7 = 5/3 and vj$ = 
7 t/ 2 . This produces a maximum value of the magnetic energy 
density of one tenth of the thermal energy density, (i.e. /? = 
10). This is still in the weak-field approximation (v_\ <C c s ). 

To follow the development of the instability in time, it was 
found useful to calculate the Lagrangian displacement field. 
This was computed by adding a vectorial tracer field x = (% x , 
X y , Xz ) to the code. The initial values of this field were simply 
the coordinates xq, yo and zq at t = 0 , and it was evolved ac¬ 
cording to the equation Dx/' l)/ = 0, making use of the same 
advection routines as in the rest of the code. The displacement 
field (0 can be found by subtracting the tracer field from its 
initial value. This Lagrangian displacement field was used for 
the representations shown in Fig. [2] 

In Fig. 0 \ve can see the evolution of the displacement field 
near the magnetic axis: it shows how a fluid volume that was 
initially a cylinder around the axis evolves in time under the 
influence of the instability of the magnetic field, in the p = 1 
case. The surfaces shown are surfaces of constant xi + Xy- 
Some magnetic diffusion (r) = 3x 10" ’’) has been added to re¬ 
move the higher spatial frequencies and give the picture a more 
pleasing look. Is it clear that we are looking at an m = 1 insta¬ 
bility since the plasma is displaced in the horizontal direction 
without losing its circular cross-section. 

In Fig. Q] we see the p = 1 and p = 2 runs compared; three 
surfaces of constant Xx + Xy have been plotted, two of which 
have been partly cut away. These confirm that in the p = 1 
case, the instability grows most quickly on the axis, and that 
in the p = 2 case, it grows most quickly away from the axis. 
This should be no surprise as here the local Alfven frequency 
wa (see Eq. (| 6 j) is greatest. 

It is possible to perform a Fourier decomposition of the dis¬ 
placements in the <j> dimension, enabling us to see the ampli¬ 
tudes of the different to modes separately. We can write the 
displacement £ as 

OO 

£(tu ,(j),z,t) = -A 0 (zu,z,t) + Z,t)e zm ^)(25) 

771=1 


and calculating 


Am Prn(/) = 


JET Iz=o A;„.A m 27rn7dn7dz \ 


TO max^z 




(26) 
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Fig. 2. Following the evolution of instabilities - surfaces of constant initial vj = 7t/10 = ct7o/5 at roughly regular time intervals: 
t = 6.30,9.07,10.86,12.52,14.56,16.20. 


should give us a suitable measure of the amplitudes of each 
m mode. A value of tu max = 0.5tuo is used, since this is the re¬ 
gion of interest, where the field strength is roughly proportional 
to zu p . In any case, choosing a different value doesn’t have any 
significant affect on the results. This quantity Amp m (or rather, 
its logarithm) is plotted in Fig.^Jfor the aforementioned p = 1 
run. We can see that the m = 1 mode is the only unstable one 
when p = 1. 

It is possible to do this Fourier decomposition on not just 
the displacement field but also on the velocity field, giving an 
analogous amplitude Amp^. This quantity is plotted in Fig. 
13 Whereas the amplitudes in the displacement field begin at 
zero, the amplitudes in the velocity field begin with a finite 
value, since we are giving the plasma an initial perturbation to 
the velocity field. 

At some point the other modes begin to grow also, pre¬ 
sumably when the m = 1 mode enters the non-linear phase. 
To test this hypothesis, we can consider that the smallest ver¬ 
tical wavelength present is two grid boxes, equal to a distance 


0.26, so the minimum value of 1 /n will be A m i n /27r = 0.042. 
The non-linear phase begins once the displacements become 
comparable to this. Fig. [3 shows the maximum value in the 
displacement field, i.e. the maximum displacement from equi¬ 
librium reached anywhere in the computational box - we can 
confirm that this value reaches A m i n /27r at roughly the same 
time as the to ^ 1 modes begin to grow. 

In the case of p = 2, the instability grows more quickly 
away from the axis. We have already seen this in Fig.0 This 
complicates the analysis, as the instability is strongest in the 
region where the initial field deviates from proportionality to 
vo p . This is to be expected, since the growth rate is of order 
v\/m which is proportional to vj if p = 2. However, it should 
be possible to use data from the inner region only, before it is 
contaminated by non-linear development from the outer region. 
Fig ,[7]is the equivalent of Fig.0]for the p = 2 case. We can see 
that the to = 0 and m = 2 modes are growing, as well as the 
to = 1 as in the previous case. 





ln[Amp m (t)] 
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Fig. 3. Two runs, identical except for the value of p, which is equal to 1 and 2 in the left and right hand figures respectively. The 
left hand figure is the same as the fourth frame of the sequence in Fig. [3 i.e. at t = 12.52. In the right hand figure, t = 15.29. 
Plotted are three surfaces of constant initial vo (i.e. constant \x + Xy) inside one another. It is clearly visible that the instability 
is on the axis in the case of p = 1, but away from the axis when p = 2. 



Fig. 4. Displacement amplitudes of different m modes, when 
p = 1. Adiabatic (k = p = 0) unstratified case. 


All subsequent discussion is limited to the physically more 
likely p = 1 case. 

4.2. Stratified and non-stratified runs and vertical 
wavenumber 

For the runs described above, we were after just a qualitative re¬ 
sult, as opposed to the quantitative result required from the fol¬ 
lowing runs, which had consequently to be set up in a slightly 
different way. 

In some cases we were interested in the behaviour at dif¬ 
ferent vertical wavelengths. We therefore performed a second 



Fig. 5. Velocity amplitudes of different m modes, when p = 1. 
Adiabatic unstratified case. 


Fourier decomposition of the displacement (or velocity) field, 
this time in the vertical direction. This yields the amplitude at 
each value of m and n (values of the latter correspond to wave¬ 
lengths from the height of the box down to two grid spacings): 

A m (tu, t) = ic m0 (ti7, t)+^2 3?(C m „(G7, t)e mz ). (27) 

n 

As before (see Eq. J26H . an integration of the coefficients C mn 
over w can be done out to a value n7 max , and it is this amplitude 
Amp mn (f) which is used for the following analysis. 

In finding the dependence of instability growth at differ¬ 
ent vertical wavelengths as a function of some parameter. 
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time 


(corresponding to a wavenumber of 20), but implemented the 
scheme described in Sect. 13.4.21 so that the gravity was in op¬ 
posite directions in the two halves of the box. A vertical res¬ 
olution of 64 was then used in a box of size 2tt x 47t/5. All 
of the following runs were done with one of these two setups. 
Since stratification has a wavelength-dependent effect, i.e. the 
stratification has a greater effect on longer wavelengths (see 
Eq.0, phenomena which either do not depend on wavelength 
or only effect the short wavelengths can be investigated using 
the non-stratified model with zero gravity. The first three of the 
following sections fall into this category: field strength, rota¬ 
tion and magnetic diffusion. As a check, the runs with rotation 
were also done with stratification. The stratified setup is how¬ 
ever required to look, in the last two sections, at the effect of 
gravity (obviously) and of thermal diffusion. 


Fig. 6. Maximum displacement, when p = 1. The computa¬ 
tional box is a cube of side 2n, and the shortest wavelengths 
of the instability may be just a twenty-fourth of this distance. 
Adiabatic unstratified case. 



Fig. 7. Displacement-field amplitudes of different m modes, 
when p = 2. Adiabatic unstratified case. 


(e.g. magnetic diffusion), one has two options: to measure the 
growth at several wavelengths in a single run with some certain 
value of the parameter in question, or to measure the growth at 
a single wavelength in several runs, each with a different value 
of the parameter. Two computational restrictions made the sec¬ 
ond method the most practical: firstly, it was discovered that at 
wavelengths comparable or greater than vjq grew more slowly 
than expected - an effect purely of the dimensions of the com¬ 
putation box, this meant that the maximum wavelength which 
could be examined was 7r/10, i.e. one fifth of voo, secondly, the 
code could not resolve very short wavelengths - the minimum 
wavelength we could look at reliably was eight grid spacings. 

In the cases where stratification was not necessary, there¬ 
fore, a vertical resolution of 8 was sufficient, with a box of 
dimensions 2tt x 7r/10; a horizontal resolution of 72 proved 
more than enough. In cases where it was necessary to hold to 
the N 3> oja condition, we still looked at a wavelength of 7r/10 


4.3. Dependence of growth rate on field strength 


The growth rate is expected to be proportional to the field 
strength. To see whether this is the case, five runs were done 
with Bq equal to 1.28, 0.404, 0.128, 0.0404 and 0.0128 respec¬ 
tively. These correspond to minimum values (that is, where the 
field is strongest) of /3 of 10 to 10 5 . Magnetic and thermal dif- 
fusivity were set to zero, as was gravity g - as described in Sect. 
I4.2l for this non-stratified case, a resolution of 72 x 72 x 8 was 
used. 

The expected growth rate is given by (from Eqs. 0 and 

0 ): 


<7 ~ UJ A 


Bo 

c7oa/ 47rp 


(28) 


Fig- |U shows the maximum growth rate reached in the veloc¬ 
ity field of this m = 1 , n = 20 mode in the five runs, against 
the value of wa- The velocity field was more convenient for 
this purpose than the displacement field, as the growth rate in 
the latter begins at infinity (because the perturbation is to the 
velocity field) and falls to a steady value, before falling again 
once the non-linear stage is reached; the growth rate in the ve¬ 
locity field begins at zero and rises to a steady value before 
falling again - its maximum value is therefore reached during 
the period of steady growth. It can be seen from Fig.[S]that the 
growth rate a is proportional to the Alfven frequency o>a. and 
is in fact almost equal to it. 

The growth rates have been predicted using the approxima¬ 
tion that the magnetic energy density is very much less than the 
thermal, and that the Alfven speed is very much less than the 
sound speed. In the run with the highest field strength these ra¬ 
tios are only 10 and 2.4, which may explain the slightly lower- 
than-expected growth rate here. 

At this point, we can have a look at the growth of the 
poloidal field component, since the production of a poloidal 
component from toroidal component via this instability is part 
of the motivation for this study. We can decompose the field 
into its three components (in cylindrical coordinates). The 
strong-field run (with Bq = 1.28, oja = 0.81) is used for 
this purpose, and the (log of the) energy contained in the three 
components is plotted in Fig. 0 i.e. a volume integrations of 
f??/87r, B^/ 87t and B\l 87rat vo < ©max = wo/2, as in 










log[max' 
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Fig. 8. The growth rate of the m = 1 mode is roughly equal to 
to a, which is proportional to field strength. Adiabatic unstrati¬ 
fied case. 



Fig. 9. The energy in the three components of the magnetic 
field. The solid line represents B i;> , dotted and dashed B z . 
The two poloidal components are growing exponentially on a 
time-scale comparable to the Alfven crossing time. Adiabatic 
unstratified case. 


Eq.|26] It can be seen in the graph that the energy in the poloidal 
components (/i„ and B z ) is growing exponentially, with a 
growth rate double that of the growth in the displacement field 
- owing to the B 2 dependence. 

In Fig. EI|is plotted the (log of the) energy in the B z com¬ 
ponent against the Alfven time, i.e. t/uo ^ 1 , for the five runs of 
different field strength Bq described above. It can be seen that 
the growth rate of the poloidal field is indeed proportional to 
ua. just as the growth rate in the displacement field in Fig. [8] 
Note that the growth rate in the strongest-field run is again a 
little lower than expected, for the reason described above. 



Fig. 10. The energy of the vertical component of the mag¬ 
netic field (B 2 /8n) against Alfven time, for five runs with 
different magnetic field strengths. The solid, dotted, dashed, 
dot-dashed and dot-dot-dot-dashed lines represent runs with 
B 0 = 1.28,0.404,0.128,0.0404,0.0128. The growth rate of 
the energy in the vertical component of the magnetic field is 
roughly is roughly equal to 2cua ■ Adiabatic unstratified case. 

4.4. Rotation 

We shall now look at the effect of rotation on the instability. 
An angular velocity is added (by means of a Coriolis force), 
parallel to the magnetic axis and gravity. In the case of p = 1 
we expect, looking at ©, to see the instability suppressed if 
fl A- to\. A number of non-stratified ( g = 0) runs were ex¬ 
ecuted, using, as before, a resolution of 72 x 72 x 8, with 
v = n = 0. The values 0, 0.15, 0.3, 0.45, 0.6, 0.65, 0.7, 0.75, 
0.9 and 1.2 of O were used, at a value of Bq of 1.28, which 
gives a growth rate in the absence of rotation of 0.64. In Fig. 
II 1 1 the amplitude of the m = 1 mode (as taken from the dis¬ 
placement field) is plotted for each of the runs, as a function of 
time. It can be seen that the instability is indeed suppressed if 
O is above a certain value somewhere between 0.65 and 0.70. 
Above this value a distinct oscillatory behaviour sets in, indi¬ 
cating stability. 

In order to check that this result holds in the regime TV A> 
to a, the runs were repeated using the stratified setup (see Sect, 
o and g = 9.1, giving TV = 5.8 so that both TV to a and 
TV A> ft. were fulfilled. This was found to make no difference 
at all to the suppression of the instability, there still existing a 
critical value of fl roughly equal to the non-rotating growth rate 
above which the instability was suppressed. 

4.5. Magnetic diffusion and its effect on high 
wavenumbers 

We expect to to see an upper limit on n due to magnetic dif¬ 
fusion, which is the same as saying that we expect to see a 
maximum value of 77 at which the instability is seen, at a given 
value of n. As explained in Sect. IQ we have to look either 
at the growth of several n modes at a given value of diffusion 
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Fig. 11. Amplitudes of the m = 1, n = 20 mode at differ¬ 
ent values of f2, taken from the displacement field. Adiabatic 
unstratified case. 


77, or to look at the growth of a particular n mode in runs with 
different values of 77 . It was found that the latter was easier 
to perform, hence we expect to see, for a given value of n, a 
value of 77 above which no instability is seen. For these runs, 
the non-stratified setup was used as gravity has no effect on 
this high-wavenumber limit. 

We used Bq = 1.28 as in the previous section, and the 
following values of magnetic diffusivity 77 : 0,2 x 10 ~ 4 ,4 x 
10 -4 ,6 x 10 - 4 ,1CT 3 ,1.6 x 10" 3 ,2.4 x 10” 3 ,4 x 1CT 3 ,6 x 
10 ~ 3 , 10" 2 and 1.6 x 10 ~ 2 . 

Measuring the growth of the mode at wavenumber n = 20 
as a function of 77, we expected to find a value of diffusion 
Writ ~ cr 0 /n 2 = 1.6 x 10 -3 (fromEq. ®) above which the 
instability does not grow, where op is the growth rate in the 
absence of diffusion. Fig.ll2lshows the growth rate of the ve¬ 
locity field measured at the time at which the growth rate in the 
zero diffusion (adiabatic) case is at its maximum (t = 10.4). It 
can be seen in the figure that the instability is not entirely sup¬ 
pressed at 77 = Writ, just slowed by a factor of two or so. Even 
when 77 is much higher than 77 cr it the growth rate is still above 
zero. 

4.6. Gravity and its effect on low wavenumbers 

We also expect to see a lower limit on the unstable vertical 
wavelength n determined by the value of g; 77 m ; n = N/uj\r 
(Eq. ®). 

For these runs we obviously have to use the stratified setup, 
the growth of the same mode, n = 20 , being measured as in 
the previous three sections. All of the runs had Bq = 0.102 
and 77 = k = 0 , each run then having a different value of gravi¬ 
tational acceleration g\ the expected value of g above which the 
n = 20 mode does not grow is Writ = nBo ^5T/8np = 10 /- 7 T; 
the code was run with the following multiples of this value: 
0,0.6,0.8,1,1.2,1.4 and 2. Fig. 1131 shows the amplitude of 
the 777 = 1,77 = 20 mode (calculated from the velocity field) 
for these values of g. As the figure shows, the instability stops 



Fig. 12. Growth rate of the 777 = 1, 77 = 20 mode in the velocity 
field as function of the magnetic diffusivity 77 . Unstratified case. 



Fig. 13. Effect of the stabilising stratification, under adiabatic 
conditions (k = 77 = 0). Velocity field amplitudes of the 777 = 
1 , 77 = 20 mode as a function of time for seven values of the 
acceleration of gravity g. The first half of the run with g = 
2 < 7 crit is plotted at higher time resolution. 


growing at the expected value g = Writ, within the measure¬ 
ment accuracy. 

4.7. Thermal diffusion 

If thermal diffusion is added, we expect to see a decrease in the 
minimum unstable vertical wavenumber, since the stable ther¬ 
mal buoyancy is reduced by diffusion on small length scales. 
To investigate this, thermal diffusion is added to the run in the 
previous section where g = Writ = 10/7T, i.e. the borderline 
case where gravity is just strong enough to suppress the insta¬ 
bility. The values of k used were 3 x 10“ 5 , 10 -4 , 3 x 10 -4 , 
1CT 3 , 3 x 10 -3 and 10 -2 . The results can be seen in Fig. M 
The figure also shows, for comparison, the case where g = 0 
and k = 0 . 
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time 

Fig. 14. Velocity field amplitudes of the m = 1, n = 20 mode 
as a function of time, at different values of thermal diffusion ft, 
when g = 10 /tt. For large values of ft, the stabilising effect of 
gravity is absent and growth occurs at the same rate as in the 
adiabatic unstratified case (uppermost line). 

From © and {QJ, we have 

2 N 2 g 2 

= u,lr*(l + n/Tr) a 1 + KnlJa (29) 

We expect, therefore, that an increase in both g and k can 
cancel each other out. It would be reasonable to expect that if 
we increase g from 0 . 8 g cr it to g CT a, we need to increase k from 
zero to ((1.0/0.8) 2 — 1 )cr/n 2 to cancel out the effect on the 
instability. This works out as k = 9 x 10~ 5 . Likewise, from 
g = 0.6<7 cr it we need k = 2.7 x 10 -4 . This means that in 
runs with g = g c r j t , using k = 9 x 10 -5 and 2.7 x 10 -4 
should produce the same instability growth rate in the m = 1 , 
n = 20 mode as if we had reduced g to 0.8g cr it and 0.6g cr it 
respectively, keeping n = 0. As Fig.ll4lshows. these numbers 
are roughly in agreement with the numerical results. 

5. Discussion 

We set out in this study to verify, with numerical methods, 
the results of the existing analytical work on the instabilities 
of toroidal fields in stars, and to check that these results were 
relevant, e.g. that a predicted instability is not drowned out by 
stronger unstable modes that might have been missed, for ex¬ 
ample due to simplifying assumptions. The results have been 
largely positive. 

Previous analytic work makes predictions of the depen¬ 
dence of instability conditions and growth rates on the param¬ 
eters of an azimuthal magnetic field. These are the dependence 
of the field strength on distance vj from the axis (the index 
p = d In B/d In w), the stability of the stratification (measured 
by the buoyancy frequency N, or equivalently the acceleration 
of gravity g ), the rotation rate O, the effects of magnetic diffu¬ 
sion (diffusivity 77 ), and the reduction of buoyancy by thermal 
diffusion on small scales (diffusivity ft). 

First we checked the analytic prediction of which azimuthal 
orders m should be unstable, depending on the value of p (Eq. 


©). We have confirmed this prediction (see Sect. 14. 1> for the 
cases p = 1 and p = 2, the former being considered the most 
important as it is this magnetic field which could plausibly be 
produced by the winding-up of a seed field by differential rota¬ 
tion. The dominant mode has m = 1, a ‘kink’ instability. 

The theory predicts that in the adiabatic, unstratified case 
(ft = 77 = N = 0) there is no threshold for instability, and 
that the growth rate a is of the order of the angular frequency 
wa = r/v a of an Alfven wave travelling around the star on 
an azimuthal field line. The field strength and hence u>a varies 
through the computational volume, but since the instability is a 
local one, the growth should be dominated by the largest value 
of u>a in the volume, after an initial transient. The numerical re¬ 
sults reproduce this well. In the best studied case, for example, 
a value cr/uj a = 0.92 was measured. 

In the adiabatic case the prediction (Pitts & Tayler 1986) 
is that the instability is suppressed when the rotation rate ex¬ 
ceeds the Alfven frequency uja ■ This was also verified, give 
or take a few percent at the most. The adiabatic case is some¬ 
what singular with respect to the effect of rotation, however. 
Theory predicts that in the presence of strong thermal diffusion 
(ft/77 — > 00) the threshold for instability disappears again, and 
that the growth rate is then of the order a = We have 

not tested this dependence, since the calculations required for 
this limiting case are computationally rather demanding. 

The effect of gravity and both magnetic and thermal diffu¬ 
sion were investigated. The effect of gravity was as expected - 
above a certain value of g the initial equilibrium is stable at a 
given vertical wavenumber n. This translates into a minimum 
unstable wavenumber which increases with increasing g. 

The effect of magnetic diffusion on the shorter wavelengths 
was not exactly the same as that expected. It was found that at 
a given wavenumber n, a value of 77 of the order suggested 
by Acheson (1978) and Spruit (1999) did not kill the insta¬ 
bility entirely, rather it reduced its growth rate by about half. 
An increase in 77 beyond this reduced the growth rate still fur¬ 
ther, but not to zero. Therefore it seems that magnetic diffusion 
alone cannot suppress the instability. It is possible that the dis¬ 
crepancy arises because the case 77 — > 00 may have a singular 
limit. The theoretical prediction was made for the case k = 0, 
but the relevant parameter is actually K/77. By analogy with 
other double-diffusive systems, instabilities must exist also in 
the case ft < 77 which do not appear when ft = 0 is set from the 
beginning. Since this case k < 77 is not of much astrophysical 
relevance, we have not pursued this further. 

Finally, the effect of thermal diffusion was tested. In the 
runs executed, it was expected that values of k of the order 
of 9 x 10 -5 and 2.7 x 10 -4 would be needed to make a run 
with g = < 7 cr it behave like the runs with g = 0 . 8 < 7 cr it and g = 
0.6(? C rit respectively. This appears to be correct, given that these 
were only order-of-magnitude approximations. 

The non-linear effect of the instability in all the above cases 
was found to be rather simple. The net effect is similar to that 
of an enhanced magnetic diffusivity: the field configuration 
spreads horizontally, while the mean azimuthal magnetic flux 
decreases due to effective reconnection across the magnetic 
axis. Toroidal fields in stars are therefore predicted to decay 
quickly by Tayler instability, once conditions for instability are 
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satisfied, unless regenerated by differential rotation. Work con¬ 
tinues on the non-linear evolution of this instability, and some 
first results can be found in Braithwaite (2005). 
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